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SUMMARY 


I A simplified two-dimensional steady- state photochemical model of the 

I atmosphere, which includes the effects of atmosphere radiation, has been 
developed to investigate the effect of injection of NO x and HO x on the 
chemical and thermal structure of the stratosphere. 

The model has been developed in two parts. The first part is the de- 
velopment of the natural atmosphere, and is characterized by the coupled 
solution of the energy and the species conservation equations in the strato- 
sphere with the temperature of the natural atmosphere prescribed. The 
results of this calculation are a coupled solution for the meridional circula- 
tion and the chemical structure of the natural atmosphere. Included in this 
model is the development of an active water vapor model that accounts for 
convection, diffusion, rainout/washout and chemistry. The model of the 
natural atmosphere has been run for the fall season, and the water vapor 
distributions are in good agreement with available data. 

The second part consists of the development of a model of the perturbed 
atmosphere and is characterized by the solution of the coupled energy and 
species conservation equations. In this model, the meridional circulation is 
prescribed from solution of the natural atmosphere and the energy equation 
is solved for the temperature perturbation resulting from the injection of 
NO and HO pollutants. The model has been used to study the effect on the 
thermal and chemical structure of the atmosphere of two types of pollution 
cases: (a) injection of NO x and HO x from a hypothetical fleet of supersonic 

and subsonic aircraft and (b) injection of HO x from a hypothetical fleet of 
liquid-fueled hydrogen aircraft. The results are discussed with regard to 
stratospheric perturbations in ozone, water vapor and temperature. 


INTRODUCTION 


The effect of aircraft pollutants on the state of the atmosphere, espe- 
cially the ozone layer, has been the subject of intense investigation for a 
number of years. In order to study this problem, photochemical models of 
the atmosphere have been developed together with making measurements 
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of various stratospheric properties. Both of these efforts have led to a 
better understanding of various stratospheric phenomena together with 
providing numerical atmospheric models that can make state-of-the-art 
estimates of the potential effect of high flying aircraft on the ozone layer. 

Atmospheric models have progressed from the most simple box model 
to more sophisticated one -dimensional and multidimensional models. Due to 
the complexity of the physical problems, these models have progressed in 
steps, developing better phenomenology for describing various atmospheric 
phenomena at each step. In all models, various aspects of the overall prob- 
lem have been either neglected or modeled in a very simple manner. In 
general, the most important phenomena have been included in these models 
and, as a result of continued development some of the initial approximations 
have been relaxed. The work reported herein attempts to relax some of 
the constraints of previous models, specifically the omission of active water 
vapor transport and the uncoupling of the stratospheric temperature distri- 
bution from changes in atmospheric species concentration. These tempera- 
ture and species concentration changes can result from the introduction of 
pollutants in the stratosphere and are of interest because of the feedback 
influence these temperature and species concentration perturbations can 
have on chemical rate processes, mean atmospheric transport and climate. 

This coupling was initially investigated by Ramanathan (ref. 1) using a 
one-dimensional radiative -convective (radiative equilibrium) model. A 
two-dimensional coupled circulation model of the stratosphere, incorporating 
the effects of radiation, photochemistry and seasonal transport processes, 
was developed by Rao (ref. 2). This model was used to study the natural 
atmosphere and to investigate ozone and temperature perturbations due to 
NO x emissions from a fleet of supersonic aircraft flying at 20 km and 45°N 
in the summer hemisphere (ref. 3). However, active water vapor was not 
included in either of these studies, and only a very simple radiation model 
was used by Rao. 

Active water vapor transport is of interest for many reasons and has 
been investigated only recently with limited success. In order to adequately 
calculate the water vapor distribution in the natural atmosphere, the dyna- 
mic interchange between the troposphere and stratosphere, together with the 
rainout /washout mechanisms in the troposphere, must be considered. In 
the initial stage of development, multidimensional models were not able to 
handle active water vapor transport; thus, the problem was circumvented by 
specifying the water vapor distributions using available measurements. 

This, however, did not allow for an active means for treating the HO x pollu- 
tion problem. This pollution problem has become much more important in 
recent years in determining the effect of aircraft exhaust emissions on the 
.chemical and thermal structure of the atmosphere. This is a consequence 
of the importance of the species OH and the measurement of new reaction 
rates that increase the relative importance of HO with respect to NO x on 
the distribution of ozone. 



The work to be described in this report consists of the development of a 
simplified two-dimensional, coupled, radiation- species model of the atmos- 
phere. Basically, it is a modification of an existing phenomenological photo- 
chemical model of the atmosphere developed at The Aerospace Corporation 
(ref. 4). Due to the limited nature of one -dimensional models, it is desirable 
to have the same technology available in a multidimensional model in order 
to more realistically determine any associated effects of pollutants. This is 
required since the more realistic transport and strong latitudinal gradients in 
temperature, density and solar flux present in a two-dimensional model can 
influence the magnitude and distribution of the resulting atmospheric per- 
turbations. 

The development of this capability requires the solution of the energy 
equation, which is a balance between the advective and turbulent diffusive 
transport of heat and the radiative heating/cooling. In this study, the local 
radiation flux is evaluated using a comprehensive simplified model developed 
by Ramanathan (ref. 1). The advective transport of heat is evaluated by cal- 
culating the meridional wind field. Since the calculation of the meridional 
circulation from the governing continuity and momentum equations is very 
time-consuming, due to numerical time step limitations, an approximate 
method is used. This method utilizes the continuity and energy equations to 
solve for the meridional velocity components. Since the energy equation is 
also utilized to calculate temperature changes, it is apparent that both the 
temperature and the wind field cannot be computed simultaneously. The pro- 
cedure in this work is to evaluate the meridional circulation in the natural 
atmosphere and then hold it fixed in the perturbed atmosphere using the energy 
equation to compute temperature changes in this case. 

In order to make these calculations self-consistent, the same radiation 
model must be used to compute both the circulation and temperature pertur- 
bations. Before the radiation model is used, a two-dimensional evaluation of 
the validity of this model is performed. The calculation of the meridional 
circulation patterns and various aspects of using this approximate method are 
also discussed. These results are compared to a similar calculation per- 
formed by Louis (ref. 5) whose circulation patterns were originally used in 
the basic model developed previous to this work. Some inert tracer experi- 
ments performed using these resulting circulation patterns are described as 
well and compared to available data. 

A discussion of the water vapor transport model is included, as well as a 
description of the coupled species-circulation-radiation model of the natural 
atmosphere. Here, the energy, species and continuity equations were solved 
simultaneously using prescribed meridional temperature distributions 
compiled from data by Louis (ref. 5). The Aerospace atmospheric model 
which was used as a basis for the incorporation of these changes is described 
in the Appendix. Typical results of this coupled calculation of the natural 
atmosphere are also presented. 

In order to evaluate the effect of HO x and NO x on the atmospheric distri- 
bution of various species concentrations and the temperature distribution, 
the meridional circulation calculated in the natural atmosphere is held fixed. 
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and the energy equation is used to solve for the associated temperature 
changes. This procedure does not allow for an estimate of the effect of the 
species concentration change on the circulation patterns, and therefore yields 
an upper limit estimate of the temperature changes. 

Two aircraft pollution cases are investigated. The first is an injection of 
NO x and HO x from a combined fleet of subsonic and supersonic aircraft flying 
in the troposphere and stratosphere. The second case considers a hypo- 
thetical fleet of liquid hydrogen-fue led aircraft cruising at 20.6 km. The 
results of these calculations are discussed in regard to perturbations in 
ozone and stratospheric temperature. 

SYMBOLS 

time constant for rainout in troposphere, s~* 
specific heat at constant pressure, cm^/s^-°K 
saturation vapor pressure, atm 
energy source/sink distribution, °K/day 
perturbation of energy source/sink distribution, °K/day 
G(z) function defined by equation (8) 

H height of troposphere 

H (</>) function defined by equation (8) 

relative humidity 

2 

turbulent eddy diffusivity tensor, km /s 

2 

coefficient of vertical eddy diffusion in troposphere, (K ), km /s 

z z 

pressure, atm 
P o reference pressure, atm 

Qg. turbulent diffusion of heat, °K/ day 

radiative heating, K/day 

Z 2 o 

R gas constant for air, cm /s' - K 

R function defined by equation (7) 

r vertical distance from center of earth, km 
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T atmospheric temperature, °K 

v meridional velocity, km/s 

w vertical velocity, km/s 

mass fraction of species i 
z vertical distance from earth surface, km 

z z/H 

(f) latitude, degrees 

( 1 potential temperature, °K 

3 

p atmospheric density, gm/cm 

reduced pressure (p/p Q )^^ C ^ 5 
X H Va o /K o 


Subscripts: 

* surface value 

Superscripts: 

( )' perturbed value 


RADIATION MODELING 


The solution of the energy equation for the temperature perturbations or 
alternatively for the meridional circulation requires the distribution of 
radiative sources and sinks in the atmosphere. In this regard, a simplified, 
yet comprehensive and realistic radiation model of the atmosphere, developed 
by Ramanathan (ref. 1) was utilized in the present study. The model con- 
siders the stratospheric radiative transfer due to t^O, CO 2 and O^. It 
includes solar absorption by six bands of water vapor in the region 0. 9p - 
b . 3 p , by ten bands of CO 2 in the regions 2. 0 p , 2. Ip and 4. 3 p , and also 
the ultraviolet and visible absorption bands of ozone. In addition, the model 
treats long wavelength radiation by the vibration- rotation, pure rotation and 
two continuum bands of water vapor, the 9. bp band of ozone and the four 
fundamental and six hot bands of CC> 2 . Specific details of the model are dis- 
cussed by Ramanathan (ref. 1). 

At the initiation of this study, this radiation model had been included 
only in a one-dimensional model, and thus it was necessary to determine 
whether the model could predict the meridional distribution of radiative 
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heating /cooling reasonably well. Asa test, the model results were compared 
with the results of a more comprehensive model developed by Dopplick 
(ref. 6). In this context, it must be emphasized that the atmospheric dis- 
tribution of radiative parameters (i. e., temperature, cloud cover, ozone and 
H 2 O distributions) was not identical in both calculations. For this verifica- 
tion study, the temperature field was taken from the data compilations of 
Louis (ref. 5), the ozone distributions from Wilcox, Nastrom and Belmont 
(ref. 7), and the tropospheric water vapor from Manabe and Wetherald (ref. 8). 
The stratospheric water vapor mass mixing ratio was assumed to be 2. 5 ppm. 

The initial results for the meridional distribution of radiative heating/ 
cooling rates were not in good agreement with Dopplick's results. After 
careful analysis of these results, it was found that the disagreement was due 
to the method in which the solar flux was evaluated. The solar heating fl ux 
i was evaluated in the original radiation model using an average daylight value 
' of the cosine of the solar zenith angle. This model does not properly eval uate | 
the diurnally averaged solar radiation needed in a two-dimensional applica- 
tion. To correct this aspect of the model, it was necessary to integrate the 
solar radiation flux over a 24-hour period. The results obtained after the 
above modification was incorporated in the model are shown in figure la and 
compared with those of Dopplick (fig. lb) for the fall season (where the fall sea- 
son refers to the northern hemisphere). Dopplick's calculations were available 
only below 32 km, whereas the present calculations were performed for the 
region 15-50 km. The present calculations were limited to altitudes above 
15 km, since the present radiation model is not valid in the troposphere. 
Comparison of figures la and lb shows good agreement between the two 
methods in the region 15-32 km. Figures 2 through 4 show the results of the 
net radiative heating rates for the winter, spring and summer seasons, 
respectively, which, although not shown, are also in good agreement with 
Dopplick’s results. 

Thus, the radiation model developed by Ramanathan provides an adequate 
description of the stratospheric radiative sources and sinks in the natural 
atmosphere, provided the diurnally averaged solar flux is evaluated properly. 
Because of its simplicity, it is well suited for use in a photochemical model 
of the atmosphere. In this regard, it can be used to evaluate the required 
radiative contribution in the energy equation for the solution of either the 
meridional circulation or the temperature perturbations. 


MERIDIONAL CIRCULATION 


The computation of the meridional circulation by solving the primitive 
equations can be very time-consuming, primarily due to the requirement of a 
small time step in order to insure that a stable numerical solution is ob- 
tained. For this reason, approximate methods have been developed. In 
general, these methods reduce the number of partial differential equations to 
two. These two equations are solved for the horizontal (v) and vertical (w) 
components of the meridional circulation, with the remaining variables 
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such as density and temperature prescribed. In most cases, one of the 
equations is the continuity equation, whereas the second is either the energy 
equation [Murgatroyd, et al. (ref. 9) and Louis (ref. 5)] or the zonal 
momentum equation [Vincent (ref. .10)]. The two-dimensional Aerospace 
Corporation model used as a basis for the present investigation originally 
employed the meridional circulation patterns computed by Louis (ref. 5). 
However, because of the differences in the radiation and turbulent diffusion 
models used by Louis and that used in the present work, the meridional 
circulation patterns must be recalculated. This is necessary in order to 
obtain a consistent evaluation of the energy equation that will allow a deter- 
mination of temperature perturbations resulting from the introduction of 
pollutants into the stratosphere. 

The method used in the present calculation follows along the lines of 
Louis (ref. 5). The major differences in the computation of the meridional 
circulation are the methods in which the radiative heating and eddy flux of 
sensible heat are modeled. The following paragraphs outline the basic equa- 
tions and the method used in computing the seasonal circulation. 


The basic equations are the energy equation 


dB 

0t 


, v 90 

r 50 


+w Ss = 



( 1 ) 


and the continuity equation 


dP + , J_ _ 

dt r cos 0 


<L 

60 


(p V C O S 0 


+ A( PW) 


=' 0 
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R /c 

where B is the potential temperature, 7 T = (p/p ) p is the reduced pres- 
sure, p is the density, v, w are the meridional and vertical velocity com- 
ponents, respectively, and T^Bn. The terms Qg ar *d Qp are the diabatic 
heating terms due to turbulent diffusion and radiation, respectively. The 
heating due to radiation is computed using the model developed by Ramanathan 
(ref. 1) a.s described in the previous section. The heating due to the turbulent 
eddy flux is computed using an eddy diffusivity model [see Reed and German 
(ref. 11)] and is given by 


6 = 1 d l *00 30) , a . ( 31 ,ii/ K dB 

r 50 \ r 60/ dz V zz 6z / r 50 V 0 Z dz 


, 6 ( K Z0 dB ) , _2_ 

6z \ r 50/ r 


k 66 , K z0 dB 
7.7. dz r 60 


tan 0 


k 90 , 90 

0 z 6z r 60 
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where IC • is the given eddy diffusivity tensor. The specific components of 

used^n this computation are those developed by Luther (ref. 12) and modi- 
fied by 7 Glatt and Widhopf ’(ref’. 13). In the above mod¥l,~we have accounted' ' - 
for both the horizontal and vertical flux contributions to the turbulent heating. 
For the calculation of the meridional circulation, the potential temperature 
and density are given. Thus, equation (1) can be written as 


+ -(Jr) 


season 


( 4 ) 


where F q represents the distribution of energy sources and sinks. As 
pointed out by Louis (ref. 5), the term dp/d t in the continuity equation is at 
least two orders of magnitude smaller than the other terms and thus can be 
neglected. Therefore, equation (2) becomes 


r cos? (PV C ° S ^ + A ^ W) " ° 


Equations (4) and (5) can now be solved for v and w. Although the compu- 
tation of v and w is straightforward, the computed values of w must 
satisfy 


L 


7T/2 


pw cos <f) d <f> = 0 


( 6 ) 


since v = 0 at both poles. 

Since the method is only approximate, equation (6) generally will not be 
satisfied. Thus, an iterative procedure is required which allows the func- 
tion F 0 to be perturbed in order to satisfy equation (6). The sequential 
steps of the solution procedure are: 

(1) Solve for w from equation (4). 

(2) Compute v from equation (5). 

(3) Evaluate the left-hand side of equation (6). 

(4) Perturb F Q by F 1 (z,0) so that the w computed by equation (4) 
with the latest values of v satisfies equation (6). 

(5) Recompute v from equation (5) using the new values of w computed 
in equation (4). 

(6) Repeat steps (4) and (5) until the values of v at consecutive itera- 
tions differ by a given €. 



Since equation (6) is an integral constraint on w, there exist many 
choices for F 1 (z ,<f>) which will satisfy this integral constraint. Thus, the 
problem of nonuniqueness of the meridional circulation can exist when this 
method is used. A sensitivity study has been performed to assess the effect 
of varying F* on the computed winds; the results will be discussed later. 


By satisfying equation (6), we require 



To simplify the integration, we assume pF'/(3#/5z) can be separated into 
the product of two functions, one a function of z and the other a function 
of 0, i.e., 

mk = G(z) H( *> ,8) 


Substituting equation (8) into the left-hand side of equation (7) and solving 
for G(z) yields 


G(z) 



cos 0 d 0 


(9) 


When equation (9) is substituted into equation (8), the following expression 
for F' (z ,<f>) is obtained. 


F' (z,0) 


l 69 
P dz 



iM 


H (<f>) c 


os 


0 d0 


( 10 ) 


The function H(0) must now be chosen. Louis (ref. 5) selected 
H (0) = cos 0. This function concentrates the perturbation around the equator 
where cos0 is unity and introduces a zero perturbation at the poles. How- 
ever, since the circulation is most important, from the standpoint of trans- 
port, at latitudes less than 50 o -60°, an attempt should be made to reduce the 
perturbation in this region. In the present analysis, H(0) is chosen to be 
unity, since for this function the perturbation is evenly distributed from pole 
to pole. The ratio of the perturbations for the function H (<f>) = cos0 (=F'jJ 
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and H{</>) = 1 (=F‘) can easily be obtained using equation (10), i.e., 

F'/F'k = 7T / ( 4 cos </>), Thus, for 4> < 45°, the perturbation using H [<p) = 1 
is less than or equal to that for the case H [(/)) = cos</>. This is desirable, 
since large perturbations are inconsistent with the method. 

This method was used to calculate the meridional circulation for each of 
the four seasons. The required inputs to the model are Q,p and the concen- 
trations of Oo, C0 2 and H^O. The seasonal distributions of temperature 
[Louis (ref. 5)], which are used in the basic two-dimensional model, were 
utilized in this calculation. The pressure was obtained from the hydrostatic 
equilibrium relationship and the density calculated from the perfect gas 
equation of state. Using the pressure and temperature, the potential tem- 
perature can be obtained. The mass mixing ratio of CO? is assumed to be 
uniform with a value of 488 ppm [Ramanathan (ref. 1)], while the monthly 
ozone distributions are taken from Wilcox, Nastrom and Belmont (ref. 7). 

The water vapor distribution is based on a relative humidity formulation 
taken from Manabe and Wetherald (ref. 8). This water vapor model was 
developed after careful analysis of data in the troposphere demonstrated that 
the meridional distribution of relative humidity was nearly similar for the 
different seasons and, thus, one could use the data to approximate an average 
relative humidity distribution. Manabe and Wetherald found the vertical 
distribution of relative humidity could be represented by 

(p/p# - 0. 02 

h = h * 1 - 0. 02 


(H) 


where h,, is the relative humidity at the earth's surface, p the pressure 
and p# £he surface pressure. Wlien p/p# is less than 0.02, the relative 
humidity becomes negative; it is therefore necessary to specify the humidity 
distribution for small values of p/p#. According to the measurements of 
Mastenbrook Cref. 14), the stratosphere is very dry and the mixing ratio is 
about 2 x 10 _fc> gm/gm of air. Thus, the mass mixing ratio is given by 

0. 622h e (T) 

Y H z O = p - h e s (T) (12) 


and 



min 


2. 0 x 10 



where e s (T) is the saturation vapor pressure at temperature T. 


Since the Ramanathan radiation model has a number of approximations 
which are only valid in the stratosphere, the model does not give an adequate 
representation for the tropospheric radiation. With a tropospheric jradiation 
model not available, the coupled continuity and energy equations 
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could only be solved in the stratosphere (i.e., above 15 km). The meridional 
circulation in the troposphere was taken from the data of Newell, et al. 

(ref. 15). This technique was also used by Louis (ref. 5). 

The meridional circulation for each of the seasons has been computed and 
the results shown in figures 5 through 8. In each case, the corresponding 
circ ulation as computed by Louis (re f. 5) is also shown. Comparison of the 
two calculations indicates similar type circulations during the winter and 
summer seasons (figs. 6 and 8) , whereas the fall a nd sp ring season s s how / 
appreciable differences (figs. 5 and 7). Since the results of these calculations 
are dependent on the modeling of the radiation, turbulent diffusion, as well as 
the method used to satisfy the integral constraint, it is therefore important to 
investigate the sensitivity of the solution to variations of modeling in these 
areas. These are shown for the fall season in figures 9 through 11 which 
depict the vertical velocity component w at 0 = 0°, 30°N and 60°N lati- 

tude. The curves labeled (a) are the results obtained using the present tech- 
nique. Curves (b) are the results obtained by Louis (ref. 5). Curves (c) are 
the results obtained using the. present model with one exception; the contri- 
bution due to eddy flux [eq. (3)] is replaced with that given by Louis (ref. 5). 
Curves (d) are the results obtained using the present model, except that in 
equation (8) H(0) = cos<£ rather than unity. Finally, curves (e).are the 
results obtained using the present model with H(<£) = cos<f> and Qj- given by 
Louis (ref. 5 ). 

Comparison.of curves (a) and (c) shows the sensitivity of the results to 
the modeling of Qjr;. The largest differences are seen to be at 30°N. Louis 
used data published by Newell, et al. (ref. 15). These data are questionable 
due to the facts that: (1) the computed fluxes are essentially the ones due to 

standing eddies only, and (2) the data in the upper stratosphere only cover 
the western half of the northern hemisphere. In the present calculations, we 
have used an eddy diffusivity model to describe the eddy flux of heat [equa- 
tion (3)]. This model accounts for both gradient and countergradient fluxes 
of heat. We should emphasize that we are not suggesting that the eddy dif- 
fusivity model predicts the fluxes of heat any more or less accurately than 
does the use of the data, but is one method of accounting for the eddy fluxes 
and is consistent with the two-dimensional diffusion coefficients used in the 
atmospheric model. 

Comparison of curves (a) [H{<£) = 1] and (d) [H(0) = cos0] shows the 
so-called "nonuniqueness " in the method of satisfying the integral constraint 
[eq. (6)]. The smallest differences occur near (f> = 30°N, since this is the 
region where the heating sources are the largest. At other locations, the 
strength of the sources and sinks are smaller and the perturbations become 
more pronounced, i.e., see (f) = 0°, 60°N. This nonuniqueness is actually a 

failure in the method and arises when the computed perturbations F' are no 
longer small when compared to F q . For the results presented, the non- 
uniqueness becomes apparent at the higher altitudes (i.e., above 25 km). 
However, the breakdown in this region is due to the poor description of the 
contribution to the heating due to the eddy flux, and thus a larger F 1 is 
required to satisfy equation (6). Comparisons of curves (b) and (e) show the 
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sensitivity of the solution to the differences between the radiation models 
used in the present calculation and that of Louis (ref. 5). As previously 
discussed, the present computation includes a radiation model developed by 
Ramanathan (ref. 1) to compute the heating due to O 3 , CC >2 and H^O, whereas 
Louis used a number of different computer codes to account for each of the 
radiating species [for details, see Louis (ref. 5)]. For the present results, 
the largest difference occurs at 0 = 60°N. In general, below 30 km, the 
magnitude of the radiation heating is less than 1 °K/day; thus there can 
exist a substantial effect if the models are predicting differences on the order 
of tenths of a degree per day w 

In general, we cannot distinctly separate out the effects of both radiation 
and eddy flux modeling due to the integral constraint imposed by equation ( 6 ). 
This is because a modification of either of the above effects will require a 
different perturbation to satisfy equation ( 6 ) and in a manner that may over- 
shadow the actual local physical effects. 

Short of solving the complete system of primitive equations for the cir- 
culation, these approximate methods are adequate and should suffice for the 
phenomenological models in which they are to be used, as long as one keeps 
in mind the sensitivity of the results to the technique in modeling the physical 
processes. Some filtering of the computed results as done by Louis (ref. 5) 
may also be necessary. 

As a further evaluation of these computed circulation patterns, inert 
tracer calculations were performed for carbon- 14 and tungsten- 185. These 
results are described in the next section. 


TRACER STUDIES 


In order to test the present meridional circulation pattern, the disper- 
sion of carbon-14 and tungsten-185 were calculated using the new meridional 
circulations computed in the previous section. Figures 12 through 14 show 
the measurements of carbon-14 distribution [ Telegadas (ref. 16); Johnston 
(ref. 17)] as a function of time for the latitudes (f> = 10°N, 30°N, and 70°N. 

Also shown are the results of the latest calculations as well as the previous 
results obtained using Louis' circulation [Widhopf, et al. (ref. 22)], At 10°N, | 

the latest circulation overpredicts the peak concentration of carbon- 14. At 1 

30°N, although the earlier profiles show the peaks occurring at 1-2 km higher I 
as well as lower concentration below the peak, these distributions approach 
the data over the next year or so. At 70°N, a similar situation occurs; i. e., 
by January 1965 the calculations adequately predict the data. 

Figure 15 shows the results of the tungsten- 185 simulation for the peak 
concentration in the equatorial region for both Louis' circulation and the 
present circulation. Along with these results are the data obtained by Friend 
(ref. 18) and Gudiksen, et al. (ref. 19). Notice that the new circulation 
adequately simulates the decay of the peak equatorial concentration of 
tungsten-185 as a function of time. 
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It must be pointed out that the turbulent diffusion coefficients used in 
these calculations were those obtained by Glatt and Widhopf (ref. 13) which 
are modifications to the original coefficients of Luther (ref. 12). This modi- 
fication involved using Louis' (ref. 5) circulation and adjusting the diffusion 
coefficients in the high latitude regions employing the carbon- 14 data as a 
guide. With these new coefficients, the tungsten- 185 simulation was then run. 
In the present calculation, the seasonal circulation was computed using these 
diffusion coefficients. Thus, any modifications to the diffusion coefficients to 
bring the present results into better agreement with the data (as was done for 
the carbon- 14) must be performed in an iterative manner since the diffusion 
coefficients couple back to the calculation of the season circulation. This 
type of study was not within the scope of the present work. However, these 
tests demonstrate the adequacy of the present transport. 

Before a coupled radiation-species-circulation model can be developed to 
study the effects of HO x pollutants on the stratospheric ozone and tempera- 
ture, a consistent model of the natural atmosphere must be developed. In 
particular, an adequate water vapor model must be used. The development 
of such a model is described in the next section. 


WATER VAPOR MODELING 


The investigation of the effects of HO x pollutants on the stratospheric 
ozone and temperature requires an active water vapor model. This model 
must account for transport by meridional winds, diffusion by turbulent eddies, 
chemical production and/or depletion, and include the effects of evaporation 
and precipitation. 

From a physical standpoint, a steady- state water vapor distribution can 
only exist if the difference between the evaporation from the surface E, and 
the precipitation P, integrated over the entire globe vanishes. Since the 
major portion of the water vapor is located in the troposphere, a tropospheric 
water vapor model must be developed. 

Before development of a two-dimensional model is initiated, one can 
make use of a simplified one-dimensional model in order to determine the 
effects of boundary conditions, transport, chemistry and rainout on the level 
of water vapor. The lowest order one -dimensional model is based on the 
assumption that under steady-state conditions in the troposphere, a balance 
exists between the divergence of the vertical eddy flux and the rainout of 
water vapor (i.e., precipitation), with chemistry being a higher order effect. 
This latter assumption is subsequently verified by actual computation. 

To elaborate on the model, consider a region of depth H in which both 
the vertical eddy diffusivity K Q and residence time A (i.e., the time 
constant for rainout) are both independent of altitude. Q Mathematically, the 
governing equation for the steady- state distribution of water vapor is given by 

d^Y 

K - A Y = 0 (13) 

O , L O 

dz 
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where Y is the mass fraction of water vapor. The boundary conditions im- 
posed on equation (13) are 


Y (0) = Y a (14a) 

»T* 

Y(H) = Y h (14b) 

where Y^ is the value of the surface mass fraction of water vapor and Y^ 
is the value of the mass fraction of water vapor in the stratosphere. The 
solution of equation (13) subject to equation (14) is given by 

Y. sinh [X (1-z)] + Y sinh [Az] 

Y = JL M (15) 

sinh A 


where z = z/H and X = H yJA /K . The first question that arises is, Does 
this model give a reasonable upward flux at surface, i. e., is the overall 
evaporation process modeled? The flux at the surface is given by 
F* = -K Q (dY/dz)^ and can be found by differentiating equation (15). The 
resulting expression is 


F . = J A K [Y. tanh A - Y„/sinh A] (16) 

O O 5:5 rl 


A computation of F* requires an estimation of the various parameters in the 
above expression. At 30°N latitude, Junge (ref. 20) estimated the average 
residence time of water vapor in the atmosphere to be about eight days, and 
the height of the troposphere is about 15 km. The value of the stratospheric 
mass fraction is about 2. 5 x 10" [see Harries (ref. 21)]. A typical value of 
K q used in_one-dimensional calculations of the atmosphere is T0 - ^ km /s 
and thus, A 6. The form of boundary condition at the surface is taken from 
Manabe and Wetherald (ref. 8) and is given by 


Y 




0. 622 h* e s (T) 
P* “ k* e s (T) 


(17) 


where e g (T) is the saturation vapor pressure, and h ;|; and p^, are the surface 
relative humidity and pressure, respectively. Typical values of Y* are on 
the order of 10"^ to 10” . Under these conditions, the surface flux can be 
approximated by 


F * = / V^o Y * 


(18) 
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When the values of A Q , K q and Y* are substituted into equation (18), the 
surface flux is found to be about 0. 2 x 10 "7 km/s. Existing data [junge 
(ref. 20)] show F* — 0.3 x 10" 7 km/s; thus the computed surface fluxes are 
considered reasonable. Figure 16 shows the results of two separate calcula- 
tions using The Aerospace Corporation's one-dimensional photochemical at- 
mospheric model along with observational data taken from Mastenbrook (ref. 
14) and Oort and Rasmussen (ref. 3 0) interpolated for 30°N latitude. The 
values of the parameters are not unreasonable and have been chosen in order 
to show the calculated distributions of water vapor that not only bracket the 
observed data, but in addition yield reasonable agreement. Curve (a) corre- 
sponds to the values: A Q = 5. 8 x 10"& day"l (two-day rainout period), K Q = 

2 x lO"^ km 2 /s, Y* = 0. 24 x 10"2 gm/gm-air; whereas curve (b) corresponds 
to A Q = 2 x 10"6 s“l (5. 8-day rainout period), K Q = 0. 66 x 10"^ km 2 /s, 

Y = 0.69 x 10" 2 gm/gm-air. These calculations give surface fluxes of about 
0.25 x 10-7 km/s. These calculations were carried out for two different 
chemistry models; namely, an original chemical model used by Widhopf, et al. 
(ref. 22) (see Table I) and that currently published by the NASA workshop (see 
Table II). Even though certain key hydroxyl reaction rates were changed, the 
results for the water vapor distributions are relatively insensitive to this 
change. The numerical results of the one-dimensional model were compared 
with equation (15) using the corresponding parameters and found to differ by 
less than two percent, thereby confirming the fact that the water vapor profile 
is mainly controlled bythe transport/rainout mechanisms in the troposphere. 


Extension of the one-dimensional model to two dimensions is generally 
straightforward. The average residence time of the water vapor as a func- 
tion of latitude is taken from Junge (ref. 20). In addition, the transport 
due to circulation as well as both gradient and countergradient turbulent 
fluxes of water vapor is included. 


Regarding the assumption of A Q being constant in the troposphere equal 
to an average value estimated from data, it should be emphasized that the 
model was not developed to predict the actual precipitation at each point in 
the troposphere (which is not known, anyway) but to simulate the total pre- 
cipitation in a column at each latitude (since the data can only show total 
precipitation that reaches the surface). In addition, since the amount of 
precipitation is proportional to the amount of water vapor present, most of 
Ifie^reTcipitatfon occurs in the first 5-6 km, and this agrees" reasonably well 
with observations of cloud cover. 


We should also point out the lack of success in attempting to develop a 
water vapor model based on a relative humidity formulation. In this type of 
model, the relative humidity is monitored; when its value reaches a certain 
predetermined level, precipitation is assumed to occur. The results of this 
type of model did not yield a dry stratosphere as observed by Mastenbrook 
(ref. 14). 


With a two-dimensional model of the water vapor available, a simplified 
coupled- radiation species and circulation model of the natural atmosphere 
can be put together. The description of the model of the natural atmosphere 
is given in the next section and includes the water- vapor distributions cal- 
culated using this model of rainout/washout. 
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NATURAL ATMOSPHERE 


In order to assess the effects of the injection of pollutants on the chemical 
and thermal structure of the atmosphere, one must have knowledge of the 
structure of the natural atmosphere. The computation of the coupled chemi- 
cal, thermal and dynamic structure of the atmosphere tends to be very expen- 
sive from a computational standpoint. To alleviate this problem, one can 
consider a simplified atmospheric model in which the perturbed thermal struc- 
ture of the atmosphere is not coupled back to the dynamics. A consistent cal- 
culation of the perturbed structure can be performed by prescribing the 
natural atmosphere by a coupled solution of the species conservation and energy 
equations with the temperature specified. The results of such a calculation 
would yield the meridional circulation as well as the chemical structure of the 
natural atmosphere. The perturbed chemical and thermal structure of the 
atmosphere can then be obtained by solving the coupled species conservation 
and energy equations, using the computed meridional circulation of the natural 
atmosphere. Although there is no coupling of the chemical and thermal struc- 
ture to the dynamics, the resultant temperature perturbation obtained can be 
considered to be an upper bound, since all the perturbed energy is converted 
into a temperature change, rather than allowing a portion of it to perturb the 
circulation. This model of the natural atmosphere is described below. 

The coupled model of the natural atmosphere is a modification of The 
Aerospace Corporation's atmospheric model. This original model was a 
steady-state seasonal phenomenological photochemical model of the atmosphere 
in which the hydrodynamic variables (mean atmospheric density, temperature, 
turbulent diffusion coefficients and mean meridional winds) are specified and 
used to solve the system of species conservation equations for the meridional 
distribution of trace species. The details of this model are described in the 
Appendix. The modification necessary to complete the model of the natural 
atmosphere was to couple the method of computing the meridional circulation 
to the solution of the species conservation equations. The resultant coupled 
model was then used to calculate the chemical and dynamic structure of the 
natural atmosphere. The water model outlined previously was also incorpo- 
rated. 

To demonstrate the feasibility of such a model, the natural atmosphere 
was simulated for the fall season (northern hemisphere). The coupled calcu- 
lation was run until a steady-state solution was achieved. Figure 17 illustrates 
the data taken from Gebhart, et al. (ref. 23), Sticksel (ref. 24) and London 
(ref. 25) for the latitudinal distribution of total ozone column for the fall sea- 
son. Also shown are the results of two calculations for the natural atmosphere. 
Curve (a) is an initial calculation carried out using the chemical system and 
associated reaction rate coefficients tabulated in Table I. This is essentially 
the chemical system used by Widhopf, Glatt and Kramer (ref. 22) to estimate 
the effect on ozone due to NO x emissions from a combined fleet of supersonic 
and subsonic aircraft projected to be operational in 1990. Curve (b) is the re- 
sult of a similar calculation in which the important reactions involving NO 3 and 
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N 2 O 5 have been included together with use of the latest updated reaction rates 
from the NASA summary (ref. 51), a s shown in Table I I. The effects of day/night 
averaging were also included in this latter calculation. The converged solu- 
tion, curve (a), shows good agreement with the data in the northern hemis- 
phere (maximum difference, about nine percent). The ozone column in the 
southern hemisphere (maximum difference, about 16 percent) is not in as good 
agreement with the available data as is the northern hemisphere calculation. 

The results given by curve (b) show a substantial increase in th_e column 
throughout the meridional plane [approximately a 50 percent increase from 
curve (a) near the northern polar region]. This result is essentially due to ■ 
the change in the reaction rate for the reaction RIO: NO + HO^-MDH +~NC >2 
from a value of 2. 3 x 10“^ to 8 x 10”^, a factor of about 35 increase. 

Since the major destruction of ozone through NO x occurs by the reaction R7, 

O 3 + NO— ►O^ + NO^, the effect of increasing RIO is to reduce the effect of 
R7. This result has also been observed in our one-dimensional model as well 
as by other investigators in both one -dimensional and two-dimensional 
models. Since the reaction rate RIO has been verified, it is possible that 
there may be other phenomena which have not been properly accounted for, 
either by having an incorrect reaction rate constant, not including pertinent 
chemical systems, or an incorrect prescription of the dynamics. However, 
to date, the latest set of reaction rates corresponding to curve (b) is the 
state-of-the-art; thus, until further refinement is made, these rates are used 
for current atmospheric model simulations. Due to the coupling of the diffu- 
sion coefficients to the computation of the circulation, it is beyond the scope 
of this work to investigate the sensitivity of the resultant ozone distribution 
to modifications in both the diffusion coefficients and chemical system. 


Figure 18 shows the vertical distribution of ozone for various latitudes 
(10°N, 40°N, 60°N and 80°N) corresponding to curves (a) and (b) in figure 17. 
The result's for solution' (a) are m agreement 'with data [Wilcox (ref. 7), and 
Hering and Borden (refs. 26 through 29)], with the peaks generally being 
about 2 km lower than the peaks in the data. The results for solution (b) show 
a large increase in the ozone level throughout the meridional plane. 

Figure 19 shows the resultant converged solution of water vapor profiles 
at 0°N, 30°N and 60°N corresponding to solution (b). Included is the cor- 
responding data of Oort and Rasmussen (ref. 30) for the tropospheric water 
vapor. The calculated results are in very good agreement with the data at 
0°N and 30°N, and slightly low at 60°N. However, these reduced values at 
60°N are most probably associated with the lower boundary value used in the 
calculation. The calculation also produces a very dry stratosphere with a 
mass fraction of 2.5-3 x 10-6 gm/gm-air. The water vapor results cor- 
responding to the chemical system in Table II are essentially the same since, 
as pointed out previously, chemistry is a higher order effect on the water 
vapor profiles. 

Figures 20 and 21 show the resultant converged solution for the meri- 
dional heating distribution and circulation pattern. Comparison of figure 21 
and figure 5 shows the computed ozone profiles give similar type streamline 
patterns in the northern hemisphere as those obtained using the prescribed 
ozone profiles by Wilcox, et al. (ref. 7). However, the circulation pattern 
in the southern hemisphere has changed somewhat. 


AT 



One important result has been in regard to the N?0 profiles. In the 
original calculation by Widhopf (ref. 4), which used the meridional cir- 
culation calculated by Louis (ref. 5), the at high altitudes was shown to 

be high when compared to data [Ehhalt, et al. (ref. 31), Farmer et al. 

(ref. 32), Murcray, et al. (ref. 33), and Schutz, et al. (ref. 34)]. These 
high values were attributed to the N 2 O being controlled by the meridional cir- 
, culation at the high latitudes, convecting it upward and then diffusing it 
meridionally. In the present calculation, the computed circulation in the 
high latitudes produces a lower upward velocity than does that of Louis and 
thus lowers the value of N 2 O. This can be seen in figure 22. 

Thus, the coupled model developed for the natural atmosphere has been 
found to produce distributions of water vapor and ozone that are in reasonable 
agreement with data. With this fact in mind, an investigation of the effects 
of NO x and HO x pollutant injection on ozone and stratospheric temperature 
can be made. 


PERTURBED ATMOSPHERE 


The development of an atmospheric model to investigate changes in both 
the chemical and thermal structures of the atmosphere, arising from the 
'injection of pollutants, involves coupling the individual spebies conser- 
vation equations to the energy equation. In this model, the circulation 
is held fixed as computed in the natural atmosphere; thus, the coupling of the 
temperature perturbations to the meridional circulation is not considered. 
Because of this lack of coupling, one can expect the calculated temperature 
perturbations to be an upper bound, since all the energy changes go into 
driving the temperature perturbation. The temperature perturbations are 
coupled back into the reaction rate coefficients in the computation of the 
chemical structure of the perturbed atmosphere. The formulation of the 
energy equation is described below. 


The energy equation expressed in terms of the potential temperature is 
given by equation (19)- 
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To simplify the calculation, equation (19) can be written in terms of the per- 
turbation of potential temperature 
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where Qj_, is given by equation (3) with 9 replaced by 9' and is - (Q„) . 
Here, the term (Qr) 0 is the resultant distribution of radiative sources ana o 
sinks obtained for the natural atmosphere, and Qp is the total radiative 
heating or cooling in the perturbed atmosphere. The term 7T 1 /tt Q (Qp) Q is, in 
general, smaller than Qk (less than one percent near the tropopause and 
about eight percent near 50 km) and thus can be neglected. Therefore, equa- 
tion (20) caii be approximated by 
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Equation (21) is solved in the domain -90°^ <f> ^ 90° and 15 ^ z 2s 50 km. 
The lower boundary of 15 km has been chosen since the radiation model 
developed by Ramanathan (ref. 1) is valid only in the stratosphere. The 
solution of the energy equation down to the surface requires the need of a 
very complex tropospheric radiation model, which is not available at the 
present time. Equation (21) requires four boundary conditions on 0'. At the 
poles we have prescribed d9' /d& = 0, whereas at the upper and lower bound- 
aries we have prescribed d9 1 /dz = 0. The use of the lower boundary condition 
will be discussed later in this section. 

The energy equation is solved using a combined leap-frog and Dufort- 
Frankel finite -difference scheme similar to that used for the solution of the 
species conservation equations. The energy equation is coupled to the species 
equations in the following manner: As the new mass fractions Y? are 
computed, they are used to evaluate the radiation term in the energy equation 
to compute 0' at time (n+l)At, The maximum time step for which the cal- 
culation was found to remain stable was about one day and was used in the 
computation of both species and temperature perturbations. With 9 ' known, 
the actual temperature T = 97T is found; thus T' = T - T Q (where T Q is the 
temperature distribution in the natural atmosphere) can be evaluated. The 
new temperature is then coupled back to the reaction rates for the computa- 
tion of the new species mass fraction. The calculation is continued until a 
steady- state solution is obtained. 

A discussion of the results of two separate investigations of the effects 
of the injection of NO x and HO x pollutants on the chemical and thermal struc- 
ture of the atmosphere follows. 
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Case 1: Perturbed Atmosphere--NO. and HO 

r xx 


Since the initiation of this study, the importance of including the effects 
of subsonic flying aircraft has become apparent. This is a result of the fact 
that the subsonic aircraft fleets were found to produce a greater amount of 
NO x than the higher flying SSTs; this potentially results in an overall increase 
of total ozone rather than a decrease [Widhopf, et al. (ref. 22)]. Table III 
shows a recent estimate (upper bound) supplied by the FAA (ref. 35) regarding 
the NO x emissions from a combined fleet of supersonic and subsonic aircraft 
projected to be operational in 1990. This table includes emissions from pres- 
ent subsonic type aircraft that nominally cruise at 9-11 km, advanced sub- 
sonic aircraft that cruise in the altitude range of 12.5-14.5 km, and super- 
sonic aircraft of the Concorde-Tupole v type that nominally cruise at 18 km. 

The variation of NO x emissions with both altitude and latitude is seen to be 
large. The major source of emissions (about 80 percent) results from the 
subsonic fleet and is deposited mainly between 9 and 13 km. The total amount 
of NO x emissions for this estimate is 2. 81 x 10^ kgm/yr. 

With the increased understanding of the HO x chemical cycle, it has become 

apparent that the effects of HO x injection are important in determining the 

chemical and thermal structure of the atmosphere. Ramanathan (ref. 1) 

studied the sensitivity of the stratospheric temperature to changes in water 

vapor using a one-dimensional radiative -equilibrium model. His results 

showed the necessity of including HO x injection for correctly predicting the 

thermal structure of the atmosphere. Two-dimensional calculations of Rao 

(ref. 3) did not include these HO effects. 

x 

The calculations to be described considered a modification of the FAA 
estimate of NO x injection rates in order to account for an increased fleet of 
supersonic aircraft. This amounted to an increase of NO x injection rates 
above 15 km by the factor 500/142. The HO x injection rates were obtained by 
multiplying the final NO x injection rates by the ratio of their emission indices 
1250/17. For the present calculation, NO x was injected as NC >2 and HO x as 
H 2 O. The initial conditions used for the calculation were obtained from the 
results of the natural atmosphere. 

The perturbed atmosphere was simulated for a period of four years after 
which a steady-state solution was attained. Figure 23 shows the latitudinal 
distribution of total ozone column change. The maximum increase in total 
ozone is three percent and occurs about 40°N latitude. This latitude cor- 
responds to that for peak NO x and HO x injection rates. Also shown, however, 
is the ozone column change above 15 km. Thus, it can be seen that more 
than two-thirds of the total column increase occurs below 15 km. The inter- 
esting result here is the production of ozone in the stratosphere. Comparison 
of the present results with those obtained by Widhopf, et al. (ref. 22), where 
NO x was the only pollutant injected, show a similar result in the troposphere; 
however, a slight reduction of ozone occurred in the stratosphere (maximum 
total ozone column change of 1.4 percent during the fall season). 
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The increase of ozone in the stratosphere in the present calculations is 
attributed to the injection of water vapor. Figure 24 shows the resultant 
distributions of H2O at (f> = 0° and 40°N for both natural and perturbed at- 
mospheres. The increase in ozone above 15 km is due to the reduction in the 
effectiveness of the NO x catalytic destruction cycle by H2O. In addition, the 
increase of the reaction rate for NO + HO2 — *-OH + NO2 can have a large effect 
on the reduced effectiveness of the NO x destruction cycle. For the pure NO x 
injection case, this catalytic cycle is 

NO z + 0( 3 P)— »NO + O z R6 

NO + 0 3 — ► NO z + O z R7 

However, the additional injection of ^O forms OH and H0 2 through reactions 

0( X D) + H z O) — ►OH + OH R32 

O s + OH— *0 2 + HO z R9 

The increased concentration of H0 2 reduces NO through the reaction 


HO z + NO— *OH + NO z RIO 

Although O3 is depleted through R9, the depletion of O3 through R7 has been 
reduced considerably due to reduction in NO through RIO. In addition, N0 2 
is reduced through the reaction 


OH + NO z + M — ►HN0 3 + M 


R12 


which in turn reduces the NO produced by the photodissociation of NO_, i. e. 


NO z + hi/— ►NO + 0( P) 


R4 


The net result is a decrease in the effectiveness of the catalytic destruction 
cycle and an increase in ozone. 

In the troposphere, the effect of the combined injection of NO x and HO 
is a production of ozone. This ozone increase is due to the increased concen- 
tration of O(^P). As pointed out by Widhopf, et al. (ref. 22) in the strict NO x 
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injection case, the increase in 0{ P) was initiated by the oxidation of 
methane by OH, which then cycled through the smog chain to produce HO 2 . 
The HO 2 reacted with NO to produce NO^ which then photodissociated to 
produce O(^P), i.e. 


HO z + NO— ►OH + N0 2 

RIO 

NO z + hi/— ►NO + 0( 3 P) 

R4 


The additional injection of HO x increases the effectiveness of smog chain on 
the production of O(^P) since H 2 0 is converted to OH by the reaction 

0( X D) + H 2 0— ►OH + OH R32 

3 

The concentration of 0( P) is also increased directly by the reaction 


OH + OH->H z O + 0( 3 P) 


R38 


3 

This increase in 0( P) increases the ozone level through the reaction 

0 2 + 0( 3 P) + M— *0 3 + M R5 

Thus, the combined effect of the injection of NO x and HO x for these cal- 
culations is to increase ozone in both the troposphere and stratosphere. It 
should be pointed out that these results are based on the new reaction rate 
measurement for RIO which has been increased by a factor of 35 over previous 
recommended value (Table I). Since the mechanisms in both stratosphere and 
troposphere are strongly dependent on the reaction HO 2 + NO — ►OH + NO?, 
and have resulted in a' substantial increase in the ozone columns in the natural 
atmosphere above observed levels, one should not draw any immediate definite 
conclusions in regard to the magnitude of the production of ozone for this NO x 
and HO x pollution case. The increased ozone levels in the natural atmosphere 
are particularly disturbing, and it is possible that other reaction rates involv- 
ing the hydroxyls are still inaccurate and could influence the present results. 
This is an important area of future investigation, and possible candidates are 
hydroxyl reactions involving production or depletion of O(^P) and those 
depleting ozone. 

Figure 25 shows the corresponding resultant steady-state meridional dis- 
tribution of temperature perturbation (T 1 ) in °K. The results exhibit two 
interesting features. First, the contours of T 1 are shown to be nearly sym- 
metrically centered around the equator in the region 40°S to 40°N between 
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15 and 25 km. Second, the maximum temperature perturbation, which is ap- 
proximately 5. 9°K, is located at the lower boundary near the equator. How- 
ever, the temperature perturbation is less than one degree above 19 km. In 
order to gain some understanding of the resultant perturbations, we must 
examine the physical mechanisms that produce these perturbations. 

Since the complete energy equation is solved in this investigation, the 
steady- state distribution of T 1 is achieved by a balance of convection, diffu- 
sion and radiation of perturbed energy. The perturbations are initially gen- 
erated through the radiation term by a perturbation in water vapor and ozone 
after which the convection and diffusion come into play by redistributing the 
perturbed energy. Since the time scales for transport are much longer than 
those for the chemistry changes, the perturbed atmosphere must be simulated 
for a number of years before a steady- state condition is attained. 

From an analysis of the steady- state solution, it was determined that the 
heating due to radiation (Qp) in the equatorial region near the lower boundary 
is produced by the 9. 6/* Og band due to the increase in ozone concentration. 

In addition, heating due to the 15/* CO 2 band occurs because of the increased 
amount of water vapor aloft which reduces the cooling by CO^ through the 
H^O overlap region. At higher altitudes, the atmosphere is cooled by the 
increased water vapor through the H 2 O emission band. The peak cooling is 
less than -0.5°K. In the equatorial region aloft, there exists a region of 
positive temperature perturbation which is due to the vertical diffusion of per- 
turbed energy from the lower regions. The nearly symmetric contours of 
T 1 in the equatorial region are a result of the closed circulation cell (see 
fig. 21) which redistributes the perturbed energy in the region 15 to 25 km. 

Since the maximum temperature perturbation occurs at the lower boundary 
and is of substantial magnitude, it is important to investigate the effect of the 
lower boundary condition on the resultant perturbations. For this study, a 
one-dimensional version of the two-dimensional model was used. In this case, 
the energy equation [eq. (21)] reduces to 
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dt dz 
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Under steady- state conditions 
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As might have been expected, equation (23) requires the total vertical flux 
of perturbed energy to be constant, i.e. 

+ F^ = constant (24) 

where F 1 = -K z (dO'/dz) and F^ = - ^(0^/77 q ) dz. We consider two differ- 
ent boundary conditions for this investigation. The first condition, 

66' /dz. = 0, is identical to the boundary condition used in the two-dimensional 
calculation, and the second is T 1 = 0 at the lower boundary. The second 
boundary condition is physically unrealistic for the present injection case 
since the major portion of the injection is below 15 km. However, some phy- 
sical insight can be gained by using such a boundary condition. 

Figure 26 shows the resultant one -dimensional distribution of tempera- 
ture perturbation for both boundary conditions, i.e., (a) h6'/hz = 0 and 
(b) T' = 0. For these calculations, we have used the injection rates at 30 u N 
latitude. Use of the boundary condition b0 ' /% z = 0 results in a vertical 
temperature distribution similar to that obtained in the two-dimensional cal- 
culation, with the peak temperature perturbation occurring at the boundary 
with a value of approximately 0. 5°K and the crossover point from positive 
to negative perturbations occurring at about 21 km. For the two-dimensional 
result (fig. 25) at 30°N latitude, the peak value of T' is approximately 1°K 
and occurs at the lower boundary with the crossover occurring at 20 km. 

Due to the large latitudinal gradient in T‘ calculated in the two-dimensional 
results, it is not unrealistic that the one-dimensional calculations predict T' 
to be about one-half that obtained with the two-dimensional calculations. 

Use of the second boundary condition T' = 0 results in a much lower 
positive temperature perturbation (T 1 « 0. 08°K) and is found to occur 1 km 
above the boundary. Above 25 km, these results show less cooling than in 
the <5.6'/<$z = 0 case; however, below 25 km, there is more cooling (or less 
heating) for the T' = 0 case. These results are consistent since the T' = 0 
case allows a finite transfer of energy downward at the lower boundary, thus 
allowing for an energy removal mechanism which did not exist for the case 
d07dz = o. 

Since the results exhibit the strong sensitivity of the temperature pertur- 
bation to the lower boundary condition, the resultant maximum temperature 
perturbation of 5. 9°K for the two-dimensional calculations would appear to 
be on the high side. In addition, as pointed out earlier, since the tempera- 
ture perturbations are not coupled back to the dynamics, one should also 
expect the magnitude of the calculated temperature perturbation to be on the 
high side. However, it is not clear at this time what the appropriate lower 
boundary condition should be. It is evident, at this'time, that the only logical 
method of alleviating the problem of the lower boundary condition is to solve 
the energy equation down to the surface and prescribe T 1 = 0 at this point. 
This method requires the use of a complex tropospheric radiation model, 
which is not available at this time, but which should be developed. 
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Case 2: Perturbed Atmosphere--H 2 0 


In recent years, there has been considerable interest in the development 
of new types of propulsion systems for high-flying aircraft. Prior to imple- 
mentation of these systems, it is important to determine whether any detri- 
mental effects to the atmosphere will arise from their use. In particular, 
the use of liquid hydrogen has been under consideration as a new type of fuel 
for these high-flying aircraft. Since the exhaust emissions from this type 
of system contain HO x , and with the recent knowledge of the increased rela- 
tive importance of HO x chemistry on the ozone layer, it is important to 
investigate this potential problem area. 

The hypothetical fleet used in this study consists of 200 liquid hydrogen- 
fueled aircraft-flying six hours per day at an altitude of 20. 6 km through the 
latitudinal corridor 35°N-55°N. The exhaust emissions were specified to be 
in the form of water vapor at a rate of 28. 5 kgm/s. These emissions were 
uniformly distributed through the corridor 35°N-55°N and between 20 and 
21 km. Starting with the natural atmosphere as the initial condition, the 
perturbed atmosphere was simulated for a period of two and one-half years 
after which a steady- state solution was essentially attained. For this calcu- 
lation, the chemical system used was that shown in Table I as modified by 
Table II. The meridional circulation was taken from the results of the natural 
atmosphere (fig. 21). 

Figure 2 7 shows the resultant meridional distribution of percent change 
in ozone column. The peak change in the northern hemisphere is about 0. 16 
percent increase, and occurs at 50°N, whereas in the southern hemisphere 
the change is about 0.28 percent increase and occurs at 70°S. In the equa- 
torial region (40°S-20°N), a 0. 12 percent decrease in ozone has occurred. 
Also shown is the column change above 15 km. Note that there is a decrease 
in ozone below 15 km except in the equatorial region where no change is 
observed. The mechanism for ozone increase in the stratosphere with regard 
to HO x injection has been discussed in the previous section. The decrease 
in ozone below 15 km (and up to 25 km in the equatorial region) is due to 
the reduction in O(^P) by H0 2 as well as by the depletion of by OH, i.e. 


H0 2 + 0( 3 P)— >OH + O z 

R 15 

0 2 + 0( 3 P) + M — *-0 3 + M 

R5 

0 3 + OH— >-0 2 + H0 2 . 

R9 


Figure 28 shows the resultant water vapor profiles at 0°N and 40°N. In the 
region of injection, the water vapor concentration has increased by about 25 
percent. The changes in total ozone column in the southern hemisphere are 
due to the water vapor being convected and diffused southward (see fig. 21). 
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Figure 29 shows the resultant meridional distribution of temperature - 
perturbation in °K. The maximum positive temperature perturbation of 
0. 7°K occurs at the lower boundary near the equator. The peak negative 
perturbation is -0.7°K and occurs near the northern polar region at about 
21 km. Note the nearly symmetrical contours centered around the equator 
below 20 km. This result is similar to that found in the NC> x and H O x 
injection, although not as wide and with a much lower zero contour. The 
cooling is initiated through the radiation term in the energy equation by the 
increased water vapor concentration after which convection and diffusion come 
into play in redistributing the energy. The heating in the equatorial region is 
due to the increased concentration of ozone as well as the increase in water 
vapor aloft which reduces the cooling of the 15/i CC >2 band through the 
transmissivity factor. Convection and diffusion then again become important 
in redistributing this energy. 

The physical aspect of the lower boundary condition with regard to the 
resulting temperature perturbations has already been discussed, and it should 
be emphasized that it also applies here. However, based on these results, 
it appears that the present hypothetical fleet of liquid hydrogen-fueled aircraft 
has a minimal effect on both the thermal and chemical structure of the atmos- 
phere. 


CONCLUSIONS 


This study was initiated to provide the ability of including active HO 
perturbations in stratospheric pollution studies and to investigate any resulting 
stratospheric temperature changes. The primary results of this investiga- 
tion are: 

(1) The Ramanathan radiation model provides a very good description of 
the meridional stratospheric heating and cooling rates, provided the solar 
radiation is calculated using a diurnal averaging technique. 

(2) The computation of the meridional circulation using the energy and 
continuity equations (given the temperature distribution) has been found to be 
sensitive to the prescription of both turbulent diffusion coefficients as well as 
the radiation model. Due to the requirement of satisfying an integral con- 
straint on the vertical velocity when using this approximate technique, the 
resultant circulations are, in general, not unique. 

(3) The resultant seasonal meridional circulation patterns computed 
using the approximate method have been used to study the dispersion of inert 
tracers (carbon-14 and tungsten- 185). The results of this study have shown 
the adequacy of these computed circulation patterns in describing the decay 
of these inert tracers. 
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(4) An active water vapor model was developed which is simple in con- \ 
cept and yields two-dimensional steady-state water vapor distributions that | 
are in agreement with available tropospheric measurements. The model also | 
yields a dry stratosphere in agreement with available data. In addition, the ! 

' calculated precipitation is in relatively good agreement with measurements, j 

being approximately 75 percent high in the region 0-30° latitude and within j 

l 25 percent of the measurements at higher latitudes. j 

| 

| (5) The inclusion of the latest chemical reaction rates has been found to j 

j produce more ozone (approximately 25 percent) in the natural atmosphere 
I than has been previously calculated and/or measured. This result indicates 
j the possibility that other chemical species are important, although presently 
i not included, or that certain reaction rates are still uncertain, in particular 
those involving the hydroxyls. 

(6) A coupled species - radiation model was developed to study active NO x 
and HO pollutant problems. It has been found herein that transport is very 
important in determining the distribution of the temperature perturbation in 
the stratosphere. Thus, one-dimensional calculations may be limited in this 

j regard. 

(7) At the initiation of this study the pollution emitted by high flying 
supersonic aircraft was to be investigated. However, subsequent work by 
Hidalgo and Crutzen (ref. 52) and Widhopf, et al. (ref. 22) indicates that the 
emissions by subsonic aircraft flying in the troposphere are important in the 
overall evaluation of the effect of aircraft emissions on atmospheric ozone. 

In this regard, the effect of NO x and HO x emissions from a modified FAA 
projected fleet of both supersonic and subsonic aircraft was investigated and 
found to produce a maximum total ozone column increase of three percent 
located at 40°N latitude. In addition, the maximum calculated temperature 
perturbation was 6°K and occurred at the lower boundary (15 km) near the 
equator. It should be emphasized that these results were obtained using the 
newly measured higher rate for the HO 2 + NO reaction which yields approxi- 
mately a 25 percent overprediction of the total ozone column- -in the natural 
atmosphere [see item (5)]. The necessity of choosing the lower boundary at 
15 km was due to the fact that the Ramanathan radiation model is only valid in 
the stratosphere. For this case, in which 80 percent of the NO x and HO x has 
been injected in the troposphere, the correct lower boundary condition for the 
temperature perturbation is uncertain. From a sensitivity study performed 
using a one -dimensional version of the two-dimensional model, the resultant 
temperature perturbations have been found to be sensitive to the lower 
boundary condition, and it would appear that the present boundary condition 
used in the two-dimensional model has overestimated the temperature pertur- 
bation. In addition, the method of computing the temperature perturbation by 
decoupling the circulation from the calculation would also lead to an over- 
prediction in the results. Thus, the 6°K perturbation is considered an over- 
estimation of the actual effect. It has become apparent that the inclusion of a 
tropospheric radiation model is a necessity in order to accurately determine 
the temperature perturbation, especially in the case of subsonic injections. 

Such a model is not available but should be developed in the near future. 
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(8) The effect of stratospheric HO x emissions by a hypothetical fleet of 
200 liquid hydrogen-fueled aircraft cruising at 20.6 km was found to be 
minimal with regard to both temperature and ozone (less than 0.25 percent 
increase in ozone and a maximum temperature perturbation of 0. 7°K). 
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APPENDIX 


Basic Atmospheric Model 


The model is a phenomenological photochemical model of the atmos- 
phere in which the hydrodynamic variables (mean atmospheric density, 
temperature, turbulent diffusion coefficients, and mean meridional winds) 
are either specified from observations or obtained indirectly from obser- 
vations as a function of time during the year and used to solve the system of 
species conservation equations for the meridional distribution of trace 
species throughout the year. The original formulation of the model, dis- 
cussed in Widhopf and Taylor (ref. 36) and Widhopf (ref. 4), is basically 
designed to examine relatively small changes in the ozone concentration as a 
function of the time of year throughout the meridional plane, since any re- 
sultant changes in the species concentration occurring as a result of the 
introduction of a pollutant are not coupled back to the atmospheric dynamics 
or temperature distributions. 

The governing species conservation equation is derived following the 
general procedure outlined by Reed and German (ref. 11) for representing 
the turbulent transport flux due to large-scale eddies. In the meridional 
plane this equation, written in terms of the mass mixing ratio, is of the form 


i^pY^ Of wY. j dpv Y.COS0 

~dt~ + ~5T~ + cos0 r<50 
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rd<}> 


dY. dY. 
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where Y^ is the mass mixing ratio P-J P of the i chemical species, p is 
the local mean atmospheric density, t is the temporal variable, r = z + R 
(where R e is the mean radius of the earth and z is the altitude measured 6 
from and normal to the earth's surface), <t is the latitude, U)- is the photo- 
chemical rate of production/depletion of the i^ 1 species, and is the local 
source/sink effect. The components of the tensor k^^j represent the diffu- 
sion coefficient in the respective directions arising from large-scale eddy 
motions, while v and w are the components of the mean circulation in the 
meridional and vertical directions, respectively. This equation is solved 
for each of the trace species considered. 
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Chemical Model 


The chemical system considered in this investigation includes the fol- 
lowing species: O^D), 0( 3 P), O z , O 3 , NO, N 2 0, N0 2 , OH, H z O, HO z , 
H^0 2 , HNO 3 , N, H, N 2 , CO, CH^, NO 3 and N 2 05 « Smog-type reactions 
initiated by the oxidation of methane by OH, which may be important in the 
lower regions of the atmosphere, are also included [Crutzen (ref. 37)]. 

These reactions also include the species CH 3 , CHO, CH 2 0, CH 3 O, CH 30 2 
and CH 3 O 2 H. The specific reaction system and the associated reaction- rate 
coefficients used in this investigation are tabulated in Tables I and II. This 
chemical system is based essentially on that recommended by the NASA-CFM 
Report. 

Computation of the absorption of solar radiation is an integral step in 
determining the chemical structure of the atmosphere, since many of the 
important reactions in the atmosphere are photochemical processes. Using 
the solar flux data compiled by Ackerman (ref. 38), the diurnally averaged 
local photodissociation rates are calculated at every location in the at- 
mosphere at each time step by a technique developed by Kramer and Widhopf 
(ref. 39). The time variation of the solar zenith angle with latitude and solar 
declination is included in the determination of J^. The absorption cross 
sections utilized to compute J. for the various species are outlined in 
Widhopf (ref. 4). 


The effect of chlorine in the atmosphere was not included in this investi- 
gation. From NAS studies (ref. 40), -the inclusion of chlorine should attenu- 
ate the effect of NO x destruction of ozone due to the formation of CjtONC^. 


Boundary Conditions 

The computational domain considered in this investigation extends from 
the north to the south pole, with a ten-degree meridional resolution; and from 
the surface to 50 kilometers, with a vertical resolution of 4z = 1 km from the 
surface to 35 kilometers, and 4z =2.5 km up to the upper boundary. At the 
polar regions, a zero latitudinal flux is assumed. 

A fixed ozone concentration [6(10)^ molecules/cc] was imposed at the 
lower boundary, as interpreted from the meridional distributions compiled by 
Dutsch (ref. 41) and Hering and Borden (refs. 26-29) [as summarized in the 
data compilation of Wu (ref. 42)]. The concentration of N 2 0 at the lower 
boundary was prescribed as an average value (0.481 ppmm) interpreted from 
the tropospheric measurements of Schutz, et al. (ref. 34), and Goldman, 
et al. (ref. 43). The latitudinal variation of the mass mixing ratio of CO at 
the surface was interpreted from the measurements of Seiler (ref. 44). The 
mass mixing ratio of CH^ (0. 75 ppmm) at the lower boundary was specified 
from the measurements of Ehhalt, et al. (ref. 31). Injection of NO and N0 2 
resulting from the anthrophotogenic activities were specified at the lower 
boundary as interpreted from the estimates of Robbins and Robinson (ref. 45). 
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The species 0( 3 P), 0(*D), OH, HO 2 , N and H were taken to be in photo- 
chemical equilibrium at the lower boundary because of their relatively short 
lifetimes, while HNO 3 , NC> 2 , NO and H 2 O 2 were removed from the tropo- 
sphere by simulating atmospheric rainout/washout. HNO 3 was removed at 
the average rates defined by Junge (ref. 46 ), while N0 2 , NO and H 2 0 2 were 
assumed to be removed at one-tenth this rate. 

The species 0( 3 P), O^D), O 3 , OH, H0 2 , H 2 0 2> N and H were assumed 
to be in photochemical equilibrium at the upper boundary, while the mass 
mixing ratios of N0 2 , N 2 0, CH^, CO and HN0 3 were analytically continued 
to the upper boundary by a second-order extrapolation in space and time de- 
scribed by Widhopf (ref. 4' ) and Widhopf and Taylor (ref. 36). This extra- 
polation allows the use of centered spatial differencing at this boundary, while 
also eliminating the necessity of specifying a boundary condition for these 
species at this location. It is an accurate and stable method of evaluating 
conditions at computational boundarie s [Victoria and Widhopf (ref. 47/)] when 
the physical mechanisms interior to the computational domain govern the 
boundary value. This is the case for N 2 O, N0 2 , CH^, CO and HNO^ which 
are being transported up into the higher regions of the stratosphere. 


Transport Data 

The meridional distributions of both the mean density and temperature 
were specified using the data obtained from ten years of observations which 
were analyzed and compiled by Louis (ref. 5). These averaged data are 
specified from the surface to 68 km for the entire meridional plane and for 
each of the four seasons. 

Luther (ref. 12) has analyzed the heat transfer, temperature and wind 
variance data of Oort and Rasmussen (ref. 30) using the procedure outlined by 
Reed and German (ref. 11) for defining the components of the anisotropic 
turbulent eddy diffusivity tensor. The three components (k^, k^z and k zz ) 
are specified for the northern hemisphere from the surface to 60 km. 

Values for the components of the diffusivity tensor in regions where obser- 
vational data were not available wer e obtai n ed by Luther by ext rapola tion, 
using the results of Hunten~(ref. 48) and Newell, et al, (ref. 49). 

These coefficients are specified for each month and were initially - used to 
parameterize the components of the turbulent diffusivity tensor. The values 
for the southern hemisphere were obtained by reflecting the northern hemi- 
spheric values, shifted by six months, and applying them appropriately in 
the southern hemisphere. However, in testing these transport coefficients 
against the dispersion of inert tracers in the atmosphere, they were found not 
to be totally adequate [Widhopf (ref. 4)] and were improved by numerical 
experimentation. The results are described by Widhopf, et al. (ref. 22). 
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Numerical Scheme 


An accurate (second order in space and time) time -dependent numerical 
finite difference scheme developed by Widhopf and Victoria (ref. 50), which 
is explicit in space and implicit in time and efficiently overcomes the 11 stiff" 
nature of the chemical system, is used to solve the governing individual 
species conservation equations. Details of the scheme as applied to this 
problem are discussed in Widhopf (ref. 4) and Widhopf and Taylor (ref. 36). 
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TABLE II. ADDITIONAL REACTIONS AND 
NEW NASA RATES 
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TAB LE III. - 1990 WORLDWIDE AIRCRAFT NO x EMISSIONS 
HIGH ESTIMATES, kg/yr 
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Oliver, R. C., et al. (ref. 35) 
















ALTITUDE 



LATITUDE 

(a) PRESENT MODEL 



LATITUDE 


(b) LOUIS (REF. 5) 

FIG. 5. FALL SEASON INTEGRATED MASS FLUX (10 12 gm/sec) 
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ALTITUDE 









ALTITUDE 





ALTITUDE 


<f> = 30°N 


(a) PRESENT METHOD 

(b) LOUIS MODEL (REF. 5) . 

(c) PRESENT MODEL WITH Q„ 
FROM LOUIS (REF. 5) 

(d) PRESENT MODEL WITH 

H(0) = COS0 

(e) LOUIS MODEL WITH Q R 
FROM RAMANATHAN 
(REF. 1) 
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VERTICAL VELOCITY, w - cm/sec 
FIG. 10. FALL SEASON VERTICAL VELOCITY 
PROFILE AT 0 = 30°N 



ALTITUDE 
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VERTICAL VELOCITY, w - cm/sec 

FIG. 11. FALL SEASON VERTICAL VELOCITY 
PROFILE AT <t> = 60°N 
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SOLID SYMBOLS: MEASUREMENT. (TELEGADAS, REF. 16) 

OPEN SYMBOLS: DATA OF TELEGADAS INTERPRETED BY JOHNSTON. ET AL (REF. 17) 

PRESENT MODEL 

WIDHOPF, ' ET AL (REF. 22) 

30 17am ,„J AN 1965 JAN 1964 OCT 1963 JULY 1963 APRIL 1963 
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FIG. 12. CONCENTRATION OF EXCESS CARBON- 14 AT 
<f> = 10°N (10^ molecule s/cm^) 





SOLID SYMBOLS: MEASUREMENT (TELEGADAS, REF. 16) 

OPEN SYMBOLS: DATA OF TELEGADAS INTERPRETED BY JOHNSTON. ET AL (REF. 17) 

PRESENT MODEL 

WIDHOPF, ET AL (REF. 22) 
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FIG. 14. CONCENTRATION OF EXCESS CARBON- 14 AT 
4 > = 70°N (1()3 molecules/cm^) 
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MONTHS AFTER INJECTION 

FIG. 15. DECAY OF PEAK TUNGSTEN- 185 IN 
EQUATORIAL REGION (dpm/SCF) 
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FIG. 18. FALL SEASON OZONE PROFILES IN THE 
NATURAL ATMOSPHERE 





FIG. 19. FALL SEASON WATER VAPOR PROFILES IN THE 
NATURAL ATMOSPHERE 



















H z O MASS MIXING RATIO (10 -6 gm/gm-air) 

FIG. 24. WATER VAPOR PROFILES FOR HO + NO INJECTION 
DURING THE FALL SEASON X X 
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LATITUDE 


FIG. 25. TEMPERATURE PERTURBATION FOR NO + HO 

INJECTION CASE (°K) DURING THE FALL^SEASC^N 
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PERCENT OZONE COLUMN 





H z O MASS MIXING RATIO (10 gm/gm-air) 

FIG. 28. VERTICAL PROFILES OF WATER FOR tLO 
INJECTION DURING THE FALL SEASON 
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